# This code replicates Table 9
rm(list=ls())
setwd("")
library(foreign)
data <- read.dta("data.dta")

# Collapse data to village level
data <- data[c("village", "treatment_catholic" , "treatment_adventist" , "treatment_mixed" , "treatment_peruana" , "treatment_maranatha", "distance_cusco", "families", "paved_place", "tourism", "lat", "lon", "elevationfeet", "distance_cusco_straight", "school_1992" , "hospital_1992", "post_station_1992", "football_field_1992", "football_field_1992_minutes", "legislature_village_president", "regularity_village_assembly", "relation_president_village", "relation_president_missionaries", "activity_shining_path")]
data <- aggregate(data, by=list(data$village), FUN=mean, na.rm=TRUE)

# Define variables
treats <- c("treatment_catholic" , "treatment_adventist" , "treatment_mixed" , "treatment_peruana" , "treatment_maranatha")
vars <- c("distance_cusco", "distance_cusco_straight", "elevationfeet", "lon", "lat", "tourism", "families", "paved_place", "legislature_village_president", "regularity_village_assembly", "relation_president_village", "relation_president_missionaries", "school_1992" , "hospital_1992", "post_station_1992", "football_field_1992", "activity_shining_path")

# Write helper function
get_stats <- function(var, treats, data) {
  c(mean(data[, var][data[, treats] ==1 ], na.rm = T), sd(data[, var][data[, treats] ==1 ], na.rm = T))
}

# Loop through variables
btable1 <- lapply(treats, function(y) t(sapply(vars, function(x) {
  get_stats(x, y, data = data)
})))

# Convert to data frame
btable1 <- data.frame(do.call(cbind, btable1))

# Adapt percentages
perc <- c("paved_place", "tourism", "school_1992" , "hospital_1992", "post_station_1992", "football_field_1992", "activity_shining_path")
btable1[rownames(btable1) %in% perc, ] <- btable1[rownames(btable1) %in% perc, ] * 100

# Add varnames
btable1$variable <- rownames(btable1)

# Add N
n <- sapply(vars, function(x) sum(!is.na(data[, x])))
btable1 <- cbind(n, btable1)

# Change colnames
cn <- paste0(rep(treats, each = 2), rep(c("_mean", "_sd"), 2))
colnames(btable1)[c(-1, -ncol(btable1))] <- cn
btable1 <- data.frame(round(btable1[, -which(colnames(btable1) == "variable")], 1))
btable1 <- round(btable1,1)

# Output table
library(xtable)
xtable(btable1,digits=1)

